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Abstract 

MicroRNAs (miRNAs) are small RNA molecules, about 22 nucleotide long, which post-transcriptionally 
regulate their target messenger RNAs (mRNAs). They accomplish key roles in gene regulatory networks, 
ranging from signaling pathways to tissue morphogenesis, and their aberrant behavior is often associated 
with the development of various diseases. Recently it has been shown that, in analogy with the better 
understood case of small RNAs in bacteria, the way miRNAs interact with their targets can be described 
in terms of a titration mechanism characterized by threshold effects, hypersensitivity of the system near 
the threshold, and prioritized cross-talk among targets. The latter characteristic has been lately identified 
as competing endogenous RNA (ceRNA) effect to mark those indirect interactions among targets of a 
common pool of miRNAs they are in competition for. Here we analyze the equilibrium and out-of- 
equilibrium properties of a general stochastic model of M miRNAs interacting with N mRNA targets. In 
particular we are able to describe in details the peculiar equilibrium and non-equilibrium phenomena that 
the system displays around the threshold: (i) maximal cross-talk and correlation between targets, (ii) 
robustness of ceRNA effect with respect to the model's parameters and in particular to the catalyticity 
of the miRNA-mRNA interaction, and (iii) anomalous response-time to external perturbations. 

Introduction 

A recently discovered molecular mechanism [l], lately named Competing Endogenous RNA (ceRNA) 
effect [2][3] , points out the importance of indirect interactions among transcript RNAs in competition for 
the same pool of microRNAs (miRNAs). MiRNAs are small - about 22 nucleotide long - non-coding 
RNAs which post-transcriptionally interact with their targets in a sequence dependent manner. In their 
mature stage, miRNAs get included in a RNA-induced silencing complex (RISC) and, eventually, thanks 
to a 6-8 nucleotide long seed region, bind specifically the miRNA response elements (MREs) in the 3'UTR 
of their target mRNAs. Depending on the degree of complementarity of the seed region, miRNAs can 
either cleave the transcripts (large overlap with the target) or downregulate their translation (low overlap 
with the target): in either case the net effect is a reduced amount of mRNAs or proteins. MiRNAs are 
known to regulate a multitude of different processes ranging from differentiation to neural plasticity, and 
their misfunctioning is often associated with the development of diseases [4][5] . 

In a nutshell the idea behind the ceRNA effect boils down to the simple observation that, while 
interacting with a target mRNA, a single miRNA cannot act on other targets. Mature miRNAs {i.e 
miRNAs loaded in RISC) are thus the limiting factor in a system of potentially interacting target mRNAs. 
If for example gene A which shares one miRNA with gene B, is up-regulated the common miRNAs will 
tend to bind preferentially to mRNA A due to its increased concentration. Consequently, mRNA of gene 
B will be less repressed resulting in a subsequent increased concentration fTU3ll6ll7^. Other studies have 
independently provided further evidences for miRNA mediated trans-regulatory mRNA effects [8j[9]. 
Since each miRNA can have several targets, a complex indirect interaction network among different 
targets emerges, where nodes are mRNA transcripts and there is a link between two nodes if they have at 
least one miRNA in common. Then, the highest the number of common miRNAs or MREs, the strongest 
the link. Such crosstalk effect has been observed in bacteria where the role of miRNAs is played by 
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small RNAs (sRNAs) and it is due to a titrative interaction among sRNAs and targets 10 . Depending 



on the number of sRNA binding elements crosstalk among sRNA targets can then be prioritized and 
selective 10 Tl]. 

Interaction via titration mechanisms entails a threshold-like behavior between the two interacting 
molecules, where the threshold position is determined by the relative amount of them 10 12-15 . This 



means that as long as the concentration of one of these two molecules is below the threshold almost all 
of them are bound in complexes with the second ones and their free amount is very low. Increasing their 
concentration beyond the threshold results in an increased amount of free molecules, while the others 
will be in turn almost all bound in complexes. Moreover, systems of molecules interacting in a titrative 
fashion also show a hypersensitivity in proximity to the threshold to changes in the molecule production 
In particular controlled conditions it has been shown that it is right near the threshold. 
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where sensitivity is maximal, that crosstalk among sRNA targets is maximal too [10, 



Remarkably, Mukherji and co-workers 16 recently observed a threshold-like effect also in miRNA 



target expression in single cells. Moreover, in line with studies in bacteria [10[14 and with earlier works on 
protein-protein interaction 12p3 , they tested a mathematical deterministic model of molecular titration 
to describe their results and found it in good agreement with experimental observations. Such results 
strengthen the idea that behind the ceRNA effect there is a miRNA-target titration mechanism. 
Motivated by 
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and [2j|3] and by results obtained in experiments with bacteria 
paper we extend previous models to the case of a general network of M miRNAs titratively interacting 
with N target mRNAs (ceRNAs) and analyze it from a stochastic point of view. So far analytical 
predictions from models for titrative interactions did not go beyond the mean- field limit [I0 plp6) or were 



limited to the case of small circuits because of the nonlinearities involved 12 . However, (i) stochasticity 
plays a central role in gene expression mostly when numbers of molecules involved are modest [l7- 19 and 
(ii) small circuits are usually embedded in more complex networks so that induced interactions might be 
relevant. Since potential crosstalk among miRNA targets is effective right in proximity to the threshold, 
where free chemical species {i.e. not bound in complexes) are present in small numbers, it is necessary a 
stochastic analysis of the system. 

Here we show that, despite the complexity and the intrinsic non-linearity of the system, a shrewd 
use of the moment generating function approach plus a simple Gaussian approximation are enough to 
obtain analytical expressions for noise and Pearson's correlation coefficients for all the molecular species 
considered in a generic network. 

As a preliminary result we describe, at the level of the independent molecular species approximation 
{viz. mean-field), the onset of a threshold-like behavior typical of titration mechanism 10 12 - 15 , which 
has been specifically investigated in ^O] in the case of a miRNA-mediated mRNA interaction, and discuss 
the possible mechanism leading to a specificity of the interactions. 

Secondly, for the first time, we derive analytical results beyond the independent molecular species 
approximation which allows for the characterization of profiles for means, noise and Pearson's correlation 
coefficients, comparing them with numerical simulations. Interestingly, we found that in proximity to 
the threshold both noise and correlation profiles among the different molecular species (in terms of 
Fano factor and coefficient of variation) show a maximum. Titration-like interactions could thus be an 
adequate mechanism to maintain homeostasis in the system: even if the noise increases, ceRNAs or 



miRNAs fluctuate in a highly correlated manner as discussed in 21-23 



Among the different parameters characterizing miRNA-mRNA interactions, the degree of catalyticity 
- i. e. the fraction of mRNA molecules that are recycled after the interaction with their target - is among 
the most disputed yet less understood ones: 
{a 
{a 



24 25 support an almost completely catalytic interaction 



' 0), while at the opposite range 26-281 support an almost completely stoichiometric interaction 
1). Finally, intermediate values of catalyticity are indeed supported by a recent work |29j. Here we 
show that ceRNA effect is robust with respect to this parameter too. In the limiting case of a completely 
catalytic interaction {i.e. 100% of the miRNA is recycled) a threshold-behavior is still observed as an 
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intrinsically out-of-equilibrium phenomenon: the location of the threshold turns out to be a monotonously 
increasing function of time such that, at equilibrium (long-time limit), no threshold behavior is observed. 

An out-of-equilibrium characteristic of the system predicted by the model is the response time of a 
ceRNA embedded in a network after a single factor perturbation. Again, in proximity to the threshold, 
we observe peculiar trends: upon switching on or off another ceRNA in the network the response times 
show a maximum and a minimum respectively, and the qualitative profiles are independent of the number 
of ceRNAs in competition. 

Finally we conclude proposing a series of specific experiments aiming at validating both qualitatively 
and quantitatively the model's predictions. 



Results 

Definition of a network of interaction miRNAs-ceRNAs 

The network we are interested in describing is schematically depicted in Figure[T]A_, where M different free 
mature miRNAs (colored stars) can interact with N different free target mRNAs (colored pentagons). 
miRNAs and target mRNAs interact via a titration- like mechanism [16| . As a first approximation we can 
think the mRNAs as irreversibly lost due to the miRNAs actions (miRNA-target association rate much 
greater than dissociation rate) while the miRNAs can eventually be recycled. Figure [lj3 shows a cartoon 
of such mechanism in which two different DNA molecules (green rectangles) are transcribed with rates 
and ka^ to become miRNA Si and mRNA Rj respectively. Eventually Si and Rj either degrade 
(broken gray stars and pentagons) with rates and gn^ or interact binding in a complex Cij via an 
effective association rate gij. 

The effective association rate gij should be thought as a combination of association, dissociation and 
degradation rates of the miRNA- mRNA complex Cij (see SI for more details). Once in the complex the 
mRNA Rj cannot be translated or utilized anymore. The parameter a (with < a < 1) is a measure 
of the catalyticity of the miRNA, that is the ability the miRNA has to be available again once having 
interact with its target. Thus, a = 1 means that for each mRNA Rj bound in a complex Cij there is also 
one miRNA Sj sequestered (and no more able to interact with its other targets) while a = implies that 
mRNA Rj effective degradation is increased by gij but this does not have any effect on the miRNA 5*^. 



Mean field approximation: threshold behavior and cross-talk 

The onset of a threshold-like response as a consequence of a titration mechanism is a rather well known 



phenomenon 10 12 16 20 . In Figures [2|\ and |3|V, we show an example of threshold effect in the case 
M = TV = 2 as a function of different ceRNA and miRNA concentrations. Such an effect can be 
derived under the assumption that the joint probability distributions of the different molecular species 
are statistically independent, as explained in Section Materials and Methods. 

In a general network of interaction of N ceRNAs and M miRNAs, when miRNA-target interaction 
strength is high, following the derivation of Eq. [Tl] and depending on the control parameter we decide 
to tune, two distinct phases emerge: (i) if all target transcription rates are below the threshold level, 
explicitly computable in terms of all other model's parameters, all targets turn out to be bounded 
in complexes and the free molecule ( i.e. not bounded) share is very low, (ii) if at least one of the 
transcription rate - say the q-th target - is above threshold, then all other target free molecule shares are 
expressed in finite amount. As shown in Figure[2]A, the emerging scenario entails a cross-talk mechanism 
where a single mRNA target above threshold is able to drive the other common mRNA targets above 
threshold. The hypothesis of a strong ceRNA-miRNA interaction can be relaxed, and still, a smoother 



threshold-like behavior is observed 10 . 



4 



Interestingly enough we note that if, as control parameter, we decide to tune the p-th miRNA tran- 
scription rate, keeping all the remaining model's parameters fixed, a mirror-like scenario emerges (as 
displayed in Figure [3]A_) : in complete analogy with the case previously discussed, also miRNAs cross-talk 
through ceRNAs. Here again, as long as all miRNAs transcription rates are below threshold, free miRNA 
molecule shares are very low. As the first miRNA transcription rate crosses the threshold, all other miR- 
NAs show a substantial increase of their free share. In this case too there is a clear cross-talk between 
miRNAs. It is interesting to note that the threshold value predicted by the model (see section Materials 
and methods) occurs at near-cquimolar concentrations of the different chemical species. 



If a hierarchy is present for the miRNA-target interaction strengths Qij / {gRiQs ) 10 20 , for example 
accounting for different miRNA regulatory elements (MREs) for different target mRNAs, then a hierarchy 
will be also established in the other target (miRNA) signal amplification levels when the amount of target 
mRNAs (miRNAs) is moved from below to above the threshold value. Targets sharing similar MREs will 



be more co-regulated than targets sharing only few MREs 20 . The miRNA-target interplay may thus 



be selective depending on the particular affinities and binding strengths 10 11 . This leads to a complex 



regulatory network with non-trivial indirect interactions among targets in competition for the same pool 
of miRNAs. 

The network sketched in Figure [T]A_ is a crude simplification of what should be a real-case ceRNA's 
network. To make things slightly more realistic see Figure |4]A_, where two groups of ceRNAs interact 
through two distinct sets of miRNAs [20]. However, a small subset of miRNAs makes the two groups 
of ceRNAs, otherwise statistically independent, weakly interacting by cross-connecting the two sets. We 
simulated the network's dynamics using the Gillespie algorithm in two different settings: in the first one, 
we modulate over time the transcription rate of one ceRNA, starting with a value below threshold and we 
first increase the transcription of one specific ceRNA (ceRNAl) rate after 35 hours. A first observation is 
that it is enough to bring above threshold a single ceRNA, to set the whole network in its non-repressed 
state. The second observation is that ceRNA-mediated regulation can be specific, i.e. we observe a clear 
hierarchy in the response of the different ceRNAs (see Figure |4j3) : those ceRNAs sharing the largest set 
of miRNA (red pentagons) respond more then the yellow pentagon set that shares a fewer number of 
ceRNAs. A second increase in the transcription rate of ceRNAl after 70 hours makes the hierarchy in 
the responses even more clear. Interestingly, also the sets of ceRNAs (orange and blue pentagons) which 
do not share any targeting miRNA respond to the over-expression of ceRNAl (although less than the 
previous two groups), thanks to an undirected effective interaction: ceRNAl pulls up the red and yellow 
pentagon sets, the yellow pentagon set pulls up the orange, and the latter the blue pentagon set. 

In the second setting (see Figure|4]C), we analyze the mirror scenario in which miRNAlO transcription 
rate is increased. Again the hierarchical responses of the different miRNAs is clearly visible. 



Beyond mean field approximation: noise and correlation coefficients 

To get insight into molecular species correlations for the miRNA-ceRNA interaction network we then 
assume that the joint probability distribution P for the different molecular species is a multivariate 
Gaussian (see section Materials and Methods). This ansatz turns out to be useful since all moments 
of a multivariate Gaussian can be expressed as a function of the first two, i.e. in terms of means 
and covariances. We will assume that the vector X = [Xi, . . . , Xn+m) '■= (Ri, ■ ■ ■ , Rn, Si, . . . , Sm) 
is distributed according a Gaussian multivariate measure of mean fi^ := E{Xi) and covariances c,;j := 
E{XiXj) ~ E{Xi)E{Xj). Thus the generic third and fourth moments read E{XiXjXk) := Cijfj,k+CikfJ.j + 
Cjkf-Li and E{XiXjXkXi) := CijCki + CikCji + cuCjk- 

In this way we are able to obtain a closed system of equations for (Xi), (Xf) and (XiXj) (see 
Supplementary Material for a detailed analysis). This assumption is not arbitrary (the usual van Kam- 
pen's expansion method [30 shows the master equation is Gaussian except for small corrections) and 
interestingly performs better than the most widely used linear noise approximation (see Supplementary 
Materials) when compared with Gillespie's simulations (see "sT] for a nice introduction to the subject). 
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Under this approximation we then find an analytical expression for means, noise and Pearson's correlation 
coefficients. 

The threshold is characterized not only by the abrupt change of the mean quantities as a function 
of the control parameter, but also by Pearson's correlation coefficients and noise (both related to the 
covariances) which turn out to show a maximum around the threshold. For each molecular species 
we evaluated in terms of variance (J(xi) ■= \/ {^f) — {Xi}"^ the Fano factor, f^^ = o■'^^ ^^/{xi), and the 
coefficient of variation, CVx^ = a(^xi)/{xi), which are both measures of noise. While the first one tells 
how much a particular process is different from a Poisson process, the second is a dispersion index. 
Figures [2j3,C and[2j3,C show such noise profiles as a function of ceRNAl or miRNAl transcription rate. 
As it is possible to notice in Figures [2j3 and |3j3, in proximity to the threshold the joint probability 
distributions are far from being independent {fxi 3> 1 for all indexes i labelling the different chemical 
species) while a multivariate Gaussian approximation is better suited to describe the simulation results. 
In Figures [2p and|3p we plot the CV profiles. Increasing the ceRNAl (miRNAl) transcription rate we 
observe a decreasing noise profile for ceRNAs (miRNAs) and an increasing one for miRNAs (ceRNAs), 
as expected because of the increasing and decreasing amount of free ceRNAs (miRNAs) and miRNAs 
(ceRNAs) respectively. Interestingly, right close to the threshold it is possible to notice a bump in the 
CV profiles. This phenomenon, due to the variances growing faster than means, is compatible with 
the bimodal distributions experimentally observed and verified via simulations in particular controlled 



conditions in bacterial sRNA target 32 , 33 



The Pearson's correlation coefficients, px x- = i^iii^ii^ — shown in Figures 



2 3 and 



3 3. The 



profile of the curves as a function of the control parameter, with a well-defined maximum, confirms the 
system hypersensitivity near the threshold. Analogously, we can define the Pearson correlation coefficient 
between miRNAs and ceRNAs (not shown). In this case, miRNAs and ceRNAs are negatively correlated. 

It is interesting to notice that exactly where the number of interacting molecules is small and the noise 
profiles show local maxima, the statistical correlation between molecular species is maximal too. Specula- 
tively, the titration interaction mechanism provides for a tool able to maintain the network homeostasis: 
potentially interacting ceRNAs (or miRNAs) needed in the same time fluctuate together . 

Threshold effect and miRNA-target catalytic interaction 

So far we considered a titrative stoichiometric (0 < a < 1) ceRNA/miRNA interaction. However, the open 
question is if cross-talk among miRNAs or miRNA targets can be possible in case of purely catalytic-like 
interaction (that is, in case of complete miRNA recycling, or rather a = in Equation [I]) |28j . 

It is straightforward to see that, at the steady state, equations for the various {Rj) (or {Si)) decouple 



when a = (see Equation 10) 20 . As a consequence, no cross-talk is possible among ceRNAs (or 
miRNAs). We found that in the out of equilibrium phase instead, the behavior is different. 

We considered the time evolution of the system in Equation 1 of the Supplementary Material, and 
then took pictures of the system at a given time t. If t is sufficiently small with respect to the time 
the complexes need to reach the steady-state, for different values of miRNA (or ceRNA) transcription 
rate we can observe the threshold behavior of Figure [5]^. Consequently ceRNAs or miRNAs cross-talk is 
possible, and statistical correlations are maximal, as shown by the Pearson's correlation coefficient profile 
in Figure [5}3. 

The emerging picture is that of a dynamical threshold whose value at a given time t tends monotonously 
to the equilibrium one in case of a 7^ and to infinity in case of a = for large time. In the latter case 
no cross-talk is observed at equilibrium (Figure [5]C,D). 

The ceRNA effect is therefore robust also in case of catalytic miRNA-target interaction, the crucial 
point lieing in the instant of time at which we look at the system. 
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Response times 

We have already discussed the threshold effect due to titrative miRNA-target interaction and how the 
system displays strong sensitivity (maximum cross-talk) and the maximal statistical correlation. We now 
want to understand how fast the system responds to an external perturbation. In particular we want 
to compute the time needed for a particular ceRNA (say ceRNAl) to reach the equilibrium after the 
instantaneous over-expression or knock-out of a second ceRNA (ceRNA2). 



Following 34 , we consider two different settings: (i) to mimic a sudden signal which saturates 
ceRNA2 promoter at i = 0, the transcription rate fc^j^ of ceRNA2 switches from zero to a given value 
(ceRNA2oFF-i.ON), (h) to mimic the opposite condition of a sudden drop of the activating signal at t = 0, 
the transcription rate of ceRNA2 k^^ switches from its initial value to zero (ceRNA2oN^OFF)- 

Defining the response time as the time needed to reach half of the way between initial and final 
ceRNAl steady state, we evaluate the response times for both switch-on (Tqn) and switch-off (Tqff) 
{i.e. for ceRNA2oFF-fON and ceRNA2oN^OFF respectively) conditions. We integrated numerically the 
deterministic system of equations obtained with M = 1 and N = 2 (see Equation 2 in Supplementary 
Material) to calculate: (i) the time Ton such that i?i(ToN) = -Rio + {Riss ^ -Rio)/2 (where Rig and Ri^^ 
are the initial and final ceRNAl steady-state respectively), (ii) the time Tqff such that Ri{Toff) — 
Rig — {Rig — Ri^J/2. The initial conditions are i?2(0) = and i?i(0) and S{0) with their steady state 
values in absence of i?2 in the former case, and R2{0) 7^ and i?i(0) and S{0) with their steady state 
values in presence of R2 in the latter. We also considered a slightly more complex network in which more 
ceRNAs are present and we compute ceRNAl response time with N = 5, 10, 20. 

We then ask two questions: (i) how the response time of ceRNAl changes at different values of basal 
miRNA concentration, and (ii) what happens when the system is complicated by the addition of other 
competing targets. 

As displayed in Figurc[6]^,B, upon increasing miRNA transcription rate ceRNAl Tqn and Tqff show 
a maximum and a minimum respectively. Both the maximum and the minimum are located at the 
threshold, where ceRNAl initial and final equilibrium values are near (see Figure [6]C). Such response 
time trend suggests an out- of- equilibrium phase transition, for which the system experiences anomalous 
dynamical features around threshold. Let us point out that around threshold, despite the change in 
terms of number of molecules from initial and final steady state is small, as depicted in Figure |6p, Tqn 
is largely increased while Tqff is decreased. Moreover, the qualitative shape of the curve is robust with 
respect to the number of targets in competition for the same miRNA (see Figure |6]4.,B where different 
line colors correspond to a different number of ceRNAs in the interaction's network): the maximum (resp. 
the minimum) of the response time depends only mildly on the number of ceRNA competitors, whereas 
the location of the threshold at which the free molecule share of ceRNAl starts being repressed depends 
linearly on the number of competitors. Moreover, the statistical correlation between ceRNAl and ceRNA2 
seems independent from the size of the ceRNA's network: the maximum level of correlation is almost 
the same upon increasing the number of ceRNAs with only a shift to higher miRNA transcription rates 
(Figure [6jl)). Therefore ceRNAl and ceRNA2 are always very correlated, notwithstanding the dynamical 
anomalies in the response-time around threshold. 



Discussion 

In this paper we analyzed the theoretical framework for the stochastic description of a general network 
of M niiRNAs interacting with N target mRNAs via a titration mechanism. With a dexterous use of the 
moment generating function approach plus simple Gaussian approximation we showed that it is possible 
to obtain analytical expressions for means and covariances for all the interacting molecules present in the 
system. 

We have first shown how the already well understood threshold effect implied by titrative interac- 
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tion 10 12 15 entails with interesting cross-talk phenomena which, so far, have been only partially 
investigated from the experimental point of view [l}|3|[7]-[9] . In particular the issue of the mirror scenario 
- for which not only ceRNAs cross-talk through competing for the same set of miRNAs, but, symmet- 
rically the same set of miRNAs too cross-talk through the common set of ceRNA - is a straightforward 
verification of the ceRNA hypothesis which, at the best of our knowledge, has never been attempted so 
far. In practice, knowing the set of miRNAs belonging to a specific ceRNA network, one could knock- 
down (resp. over-express) a given miRNA in the network. In this case, the model predicts that the 
other miRNAs in the network, driven by the controlled miRNA knock-down (resp. over-expression), 
should increase (resp. decrease) their free molecule share. Such an effect could be directly measurable 
as a down-regulation (resp. up-regulation) of any of the miRNAs targets (either belonging to the same 
ceRNA network, or to any other secondary target). 

In addition to cross-talk and threshold phenomena, the model predicts interesting and experimentally 
measurable trends for the noise and Pearson's correlation coefficient profiles. In proximity to the thresh- 
old, where all the free molecular species involved in the system are present in small numbers, both the 
noise measures we analyzed (Fano factor and coefficient of variation) show a maximum (for the latter 
coefficient the maximum is local). These behaviors are interpretable in terms of bimodal distributions 
for each molecular species involved in the titrative mechanism [33| . Interestingly the bimodality has been 



experimentally measured in a simple sRNA-mediated circuit in Bacteria 32 , and could be potentially 
verified in our ceRNA case. 

In proximity to such threshold value, also the Pearson's correlation coefficients among ceRNAs or miR- 
NAs show a maximum, meaning that the statistical correlation among molecules deriving from different 
genes is high. That is, not only the system is hypersensitive to little changes in the control parameter, 
but also fiuctuations are highly correlated. As a matter of fact, the titration mechanism of interaction 
establishes a positive coupling among ceRNAs belonging to different genes (or among miRNAs). While 
the intensity of such correlation depends mostly on the combination the basal transcription rates of each 
particular gene (so that different genes speak each other at different intensities, but the level of corre- 
lation is established by the particular parameters), the location of the maximum is a determined by all 
the molecular species in competition. Furthermore, such statistical correlation is robust with respect to 
the number of ceRNAs involved in the system (with just a shift in the location of the threshold when 
increasing the number of ceRNAs) and also with respect to the catalyticity parameter a. When a is 
zero, meaning that all the miRNAs are recycled, it is still possible to observe the threshold effect and the 
maximum in correlations' profiles as an out-of-equilibrium characteristic of the system. Thus, the ceRNA 
effect is always present, provided that the observation's time is short enough. 

To investigate experimentally these features, quantitative fluorescence microscopy seems, for the time 



being, the most promising technique. Previous works not directly related to the ceRNA hypothesis (see 10 



for a seminal work in bacteria, and 16 in human cell lines) used two-colors fluorescent reporter systems. 
The construct typically consists of a bidirectional drug-inducible promoter driving the expression of the 
two fluorescent proteins. The 3'UTR of the fluorescent proteins can be engineered to control the binding 



sites, and so the miRNA-mRNA binding affinity for the targeting miRNAs of interest. Both in 10 



and 16 , the method was used to monitor the threshold effect in a simple sRNA/miRNA — > mRNA 
interaction. At the expenses of creating more complex constructs, an analogous technique could be 
deployed to investigate threshold, cross-talk, and noise/correlation behavior of simple ceRNA networks. 
In the most straightforward implementation one needs two reporter constructs: (i) the first construct 
consists of a bidirectional reporter system composed by the 3'UTR of ceRNAl concatenated to the 
fiuorescent gene (say green), and on the other side a miRNA binding site free 3'UTR concatenated to a 
second fluorescent gene (say yellow) to monitor the transcription activity, (ii) the second construct consists 
of a single reporter composed by the 3'UTR of ceRNA2 concatenated with a third fluorescent gene (say 
cherry) . In this way one could simultaneously monitor the activity of both ceRNAs (green, cherry) as a 
function of the transcriptional activity of ceRNAl (yellow) which would validate both qualitatively (in 
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terms of the profile predicted by the model) and possibly quantitatively (by allowing a multi-parametric 
fit of the model's kinetic constants from the experimental data) the model predictions as displayed, for 
instance, in Figure [2] 

Finally, the model shows interesting out-of-equilibrium features around threshold which could be 
experimentally testable. In particular the peculiar response time profile as a function of the distance from 
the threshold could be directly measured by means of quantitative time-lapse fluorescence microscopy \3^ 
and flow cytometry to monitor ceRNAs dynamics. To monitor the dynamics of two ceRNAs, one could 
conservatively construct a two colors fluorescent reporter system that allows for simultaneous monitoring 



of protein levels (see again 10 16]). Of course larger networks could be potentially monitored using 
multiple colors. 

Materials and Methods 
Stochastic simulations 



Stochastic simulations have been performed via implementation of Gillespie's first reaction algorithm 36 



Theoretical framework: stochastic model 

In analogy with Figure[lj3, for each gene belonging to the miRNA-target network in Figure[T|\ we consider 
the key steps of transcription, degradation and titrative interaction among transcripts. Thus, the system 
is described by M -I- TV variables (M miRNAs Si and N target mRNAs Rj transcribed from M + N 
different genes) and the probability of finding in a cell exactly R, S := S*!, . . . , Sm, Ri, ■ ■ ■ , Rn molecules 
at time t satisfies the following master equation: 

M N 

dtP = J2ksAPs^-i-P)+Y,kRAPR,-i-P)+ (1) 

i=l j=l 

M N 

+ E 9s. m + l)Ps.+i - S,P) + E gn, {{R, + 1)Pb,+i - RjP) + 

M N 

+ aY,T.9^AiS^ + + l)Ps^+iM,+i - S.RjP) + 

J=l j=l 

M N 

+ (1 - EE 9^,S^{{R, + l)Pfl,+l - RjP) , 
i=l j = l 

with P = Pxi,...,Xk:....XM+N and Pxk±i = Pxi,...,Xk±i,...,XM+N ■ In Equationfljfcs, and k^. are transcrip- 
tion rates and 35. and g^^ degradation rates for the i-th miRNA and the j—tn target mRNA respectively. 
gij is the effective association rate for miRNA Si and its target Rj. a is the catalyticity parameter de- 
scribed above. 

By defining the generating function, 

M TV 

F(z,q|t)=Enn^f"9f (2) 

S,Ri=l j = l 

where z, q := zi, . . . , zm, li, ■ ■ ■ 1 Qn i we can convert Equation [T] into the following second-order partial 
differential equation: 

dtF{z,q\t)^Hiz,q)F{z,q\t) (3) 
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where the operator 7^(2, q) is defined as: 

M N 

H(z,q) = ^fcs.(z,-l)+^fci^^.(g,-l)+ (4) 

1=1 j=i 

M N 

+ 9s, {dz, - Zid^^ ) + gn^ {dq^ - qjdg^ ) + 

M N M N 

i—1 j — 1 i—1 j — 1 

The moment generating function has the fohowing properties: 

F(z = l,q=l) = 1, (5) 

dziP\z=l,q=l = {Si) , 

dlF\,=i,^=^ = (^■>-(^.>, 

= {S,R,). (6) 

Considering higher order derivatives in Equation |3] at steady state {dtF — 0), and assuming that all 
derivatives are computed in z = l,q = 1, we find; 

fcs, — aY\^. Qi-jd'?. „.F 
(S,) = d.,F=— '""^ , (7) 



3s 

(Rj) = dq^F 



kn, - Ei=i9ij9l,q,F 



9R, 

ksA^ + dz.F) ~aY.U 9^3^91. ^F + dl^g^F) 

(Sf) = dlF + d,.F= — 

.9s. 

(1 + dq^F) ~ Efii g.Adl q^F + dl^^F) 

n — ' 



= dlF + dq^F 



kBA.P + ksA^F - T.fii 9ijdl..„q,F- a Y.t, 9udl,q„qF 



{SrRj) 

etc... . 



9ij + 9 Si + 9Rj 



The moment-generating function defined in Equation [3] is unfortunately too complicated to be computed 
analytically even at steady state, as all moments depend on higher ones and the system is not closed, as 
shown in Equation [Tj In the following we will present a series of increasingly accurate approximations 
for analyzing it. 
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Independent molecular-species approximation 

As a first step for determining analytically the behavior of the system, we will assume that the probability 
distribution P is factorized: 

M N 

P'"^(R,S) ■.= l[Pf{S.)l[Pf{Rj) (8) 

Under this assumption it turns out that the steady state solution for the Pf{Si), and Pj^{Rj) are Poisson 
distributions whose mean value can be expressed solving the following second order system of equations, 

, , fcs, -a(5'j)indEj^i3y(^j>ind .n^ 

(5j)i„d = 1 < I < M (9) 

9s, 

I n \ fcfl., - (-Rj)ind Y!lL\9i3kSi)\^& 1 ^ ■ ^ Af 

\-Kj/ind = 1 S J S JV . 

9R, 



Analytic solutions for the system of equations 10 can be easily written in the case gji- ~ gn., 9Si = 9s 
and Qij = g for all Rj and Si'. 

1 / N M r-r 

iv \ '^^j> t I, , I. 1 gags- V A 

with A = 4:ggsgRaJ2j^^kii. + {gngs + giYlfii^Si - a^jLi^flj))^ • In the more general and 
biologically relevant case of different molecules half- lives and complex affinities gij , solutions can still be 
found, but they turn out to be too complex and long to be reported here. 

Locating the threshold 

The simplest way to locate the threshold is to solve the system of equations [TO] in the limit of strong 
interaction miRNA-target (high gij) thus finding: 



If «E,=ifc«, <E,=ifc5. (11) 

otherwise 

{ otherwise 

The threshold position is determined by the relative amount of miRNAs and their targets (see Equa- 
tion 11 ). For fixed kji- and ks^ , with j = {1, ...,q — l,q + 1, N} and i — {1, M}, the threshold 
is setby ka^ and by all miRNA transcription rates ks^- Thus, as long as the q-th mRNA target tran- 
scription rate kn^ is below its threshold level fc|j = {J2fLi ^Si — ce Y^^^q kn^ )/a all targets are bound in 
complexes and their free molecule amount is very low (while miRNAs are expressed), or, in other terms, 
the threshold is located at near-equimolar concentration of the different chemical species. 

Increasing fc^^ beyond its threshold results in the expression of all the other targets (while miRNAs 
will be all bound in complexes), see Figure [2K. 
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Within the independent chemical species approximation in Equation [8] the Fano factor (noise index 
fi^X) = crp^^/(X)) for each molecular species is 1. The factorized approximation is good enough in 
showing the threshold effect, but fails in determining correlations among molecular species (see symbols, 
which are the results of Gillespie's simulations, in Figures [2]A_ andjsj'V). 



Gaussian Approximation 

The simplest approximation beyond mean-field is a Gaussian one. We denote X = (Xi, . . . jXn^hi) := 
(i?i, . . . , Rn, Si, ... , Sm)- The approximation assumes that X is distributed as a multivariate Gauss: 



P{X) 



exp 



i{x-firc-\x~fi) 



^(27r)W+Mdet(C) 



(12) 



where the covariance matrix C has coordinates Cij :— E{XiXj) — E{Xi)E{Xj), the vector fl has coor- 
dinates :— E{Xi), and the expectation value E(-) is with respect to the Gaussian measure P defined 
in Equation |12[ All moments of a Gaussian multivariate measure can be expressed in terms of fii and 
Cij . Therefore the moments derived from the generating function in Equation [7] can be expressed in 
terms of fj,i and Cy . In the Supplementary Material we describe in details the computation of the specific 
N = M = 2 case, and we compare the performance of the Gaussian approximation with the linear-noise 
approximation. 



Acknowledgments 

While completing this manuscript we learned that M. Figliuzzi, E. Marinari, and A. De Martino have 
independently studied the same problem, reporting results which are consistent with those obtained here. 

We thank Michele Caselle, Enzo Marinari, Paolo Provero, Andrea De Martino, Luca Dall'Asta, Carlo 
Baldassi, Matteo Osella, Marco Zamparo, and Matteo Figliuzzi, for interesting discussions about technical 
aspects of stochastic modeling. We are indebted with Pier Paolo Pandolfi, Yvonne Tay, Florian Karreth, 
and Riccardo Taulli for many illuminating discussions about the experimental strategies for validating 
the model, and Terence Hwa for pointing us a relevant bibliographic reference on the subject. RZ 
acknowledges support from the ERG Grant No. OPTINF 267915. 



12 




Figure 1. Representation of a generic miRNA-target interaction network. (A) Simplified picture of a 
miRNA-ceRNA interaction network. (B) For each miRNA (Si) and ceRNA (Rj) present in the network 
we consider the main steps of transcription (rates ^5. and fc^^. respectively) and degradation (rates gs^ 
and gu- respectively) plus a titrative interaction between miRNA and ceRNA. miRNA and ceRNA can 
therefore form a complex Cij with effective association rate gij . The parameter a (the catalyticity 
parameter) tells which is the probability a miRNA is recycled after having interact with one of its 
targets. 
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Figure 2. Threshold, noise and Pearson's correlation coefficient varying ceRNA transcription rate. 
(A-C) Steady state value for means, Fano factors and coefficients of variation for each free molecular 
species in a system of two miRNAs (miRNAl and miRNA2, green and orange lines respectively) 
interacting with two ceRNAs (ceRNAl and ceRNA2, blue and red lines respectively) varying the 
concentration of ceRNAl. In proximity to the threshold the system shows hypersensitivity to changes 
in the control parameter (ceRNAl transcription rate), captured by a maximum in the Fano factors 
(panel B). For the same values of ceRNAl transcription rate, the local maximum in the coefficients of 
variation (panel C) is the fingerprint of bimodal distributions in the number of molecules for each 
molecular species. (D) Pearson's coefficients between the two miRNAs (orange line) and the two 
ceRNAs (blue line). The two lines show a maximum in proximity to the ceRNAl transcriptiom rate 
threshold value, meaning that there is a region of parameters where the fluctuations in the number of 
ceRNAs or miRNAs are highly correlated. Lines are the results of Gaussian approximation while 
symbols are Gillespie's simulations. For panels B,C the line color-code is the same as in panel A. 
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Figure 3. Threshold, noise and Pearson's correlation coefficient varying miRNA transcription 
rate. (A-C) Steady state value for means, Fano factors and coefficients of variation for each free 
molecular species in a system of two miRNAs (miRNAI and miRNA2, green and orange lines 
respectively) interacting with two ceRNAs (ceRNAl and ceRNA2, blue and red lines respectively) 
varying the concentration of miRNAI. In proximity to the threshold the system shows hypersensitivity 
to changes in the control parameter (miRNAI transcription rate), captured by a maximum in the Fano 
factors (panel B). For the same values of miRNAI transcription rate, the local maximum in the 
coefficients of variation (panel C) is the fingerprint of bimodal distributions in the number of molecules 
for each molecular species. (D) Pearson's coefficients between the two miRNAs (orange line) and the 
two ceRNAs (blue line). The two lines show a maximum in proximity to the miRNAI transcriptiom 
rate threshold value, meaning that there is a region of parameters where the fluctuations in the number 
of ceRNAs or miRNAs are highly correlated. Lines are the results of Gaussian approximation while 
symbols are Gillespie's simulations. For panels B,C the line color-code is the same as in panel A. 
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Figure 4. Selectivity of miRNA and ceRNA interactions. (A) Example of a network of ten miRNAs 
interacting with ten ceRNAs in blocks. The interaction links are such that we can define two main 
blocks (block 1 and block2) of strongly interacting miRNAs-ceRNAs connected by two common 
miRNAs (miRNA 4 and 5 in block 1, miRNA 6 and 7 in block 2) and ceRNAs (ceRNA 4 and 5 in block 
1 and ceRNA 6 and 7 in block 2). Panels (B,C) show an example of dynamics of such network. Varying 
ceRNAI (panel B) or miRNAlO (panel C) transcription rate during time (every 35 hours in the 
example, but the time is arbitrary) has a differentiated effect on the other ceRNAs and miRNAs present 
in the all network. The color-code for lines in panels B and C follows the color of miRNAs and ceRNAs 
in panel A. 
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Figure 5. Threshold effect in a miRNA-target catalytic interaction. Example of a system of one 
miRNA interacting with two ceRNAs with catahticity parameter a = 0. The threshold effect is possible 
only if the system is out of equilibrium (A). Numerical integration of Equation (1) in Supplementary 
Materials leads to time evolution of each molecular species for a given set of parameters. In panels A,C 
we plot "pictures" of the evolving system at different time t (panel A t = 10'^, panel G t = 10^) as a 
function of ceRNAl transcription rate. When t is smaller than the time complexes need to reach their 
steady state a threshold effect is observed. In panels B,D we plot the corresponding Pearson's 
coefficient profiles. Where the threshold effect is present (panel A), a peak in the Pearson's coefficient is 
also observed. 
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Figure 6. Response times upon one ceRNA perturbation.lncTeasing miRNA transcription rate ceRNAl 
shows a maximum and a minimum in its response times upon switching on or off ceRNA2 transcription 
respectively (panel A and B). The maximum (minimum) is located near the threshold, where ceRNAl 
initial value (that is its values before switching on (off) ceRNA2) is near to the steady state it will reach 
upon switching on (off) ceRNA2 (panel C) but also the more sensitive to ceRNA2 variation (look at the 
maximum in the Pearson's correlation coefficient between ceRNAl and ceRNA2 in panel D). Different 
color lines correspond to different numbers of ceRNAs in competition for the same miRNA. The 
qualitative trend for response times and Pearson's correlation coefficient is robust with respect to 
increasing such number. 
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A Generalized mean-field equation with explicit complexes 

We describe the general case of N different target mRNAs interacting with M different miRNAs. The 
action of a miRNA on its target has the following characteristics: each miRNA molecule can constitute 
a complex with a target molecule and then can be eventually released. The molecular species are: 
free miRNAs {Si), free mRNAs (Rj), complexes Cij of miRNA Si with niRNA Rj. Each gene can be 
transcribed with rate k^jf. g.j, degraded with rate g{R-^Sj}- Complexes associate with rate fc^ and 
dissociate with rate k~-. Each complex eventually degrade with rate jij. A schema of such network is 
represented in Figure m The mean-field equations thus reads: 

dRi ^' 

= fcfl. - 9R.Sj + ^ {-k+jSjR^ + kljdj) 



dt 

dS 
~dt 



N 



' - ks, - gs, Sj + J2 {~k+S,Rj + k^f,, + (1 - cy)^^,C,, 



dC 

= k+R.S,^{k^+^,i)C,, (13) 
with j £ {i, . . . , Af } and i £ {1, . . . , N}. Assuming that complexes reach the equilibrium faster than the 



other molecular species, we can simplify the system 13 to the following one 

dS, 



ksi - 9S,S^ - ag.jSiRj (14) 



dRj 



kRj — QRjRj — OijSiRj 
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Figure 7. Representation of a generic miRNA-target interaction network. (A) Simplified picture of a 
miRNA-ceRNA interaction network. (B) For each miRNA [Si) and ceRNA {Rj) present in the network 
we consider the main steps of transcription (rates ^5. and kg. respectively) and degradation (rates gsi 
and respectively) plus a titrative interaction between miRNA and ceRNA. miRNA and ceRNA can 
therefore form a complex Cij with association rate fc^^and dissociation rate k~-. The complex can then 
degrade with rate 7^^ . The parameter a (the catalyticity parameter) tells which is the probability a 
miRNA is recycled after having interact with one of its targets. 

B Generalized master equation with explicit complexes 

The master equation corresponding to Equation [T3| reads: 

M N 
dtP = X! '^sAPSi-l,Rj,C,j ~ PSi,Rj,C,j) + X! kRj{PSi,Rj-lfii, - PSi,R,,C,,) + (15) 

M N 

+ X! 3Si {{St + l)Psi+i,R^,c,j ,Rj + l,C,j - P-jPSi,Rj,Ci.j) + 

i=l j=l 
M N 

,+i,i?j+i,c.j-i - S^RjPsi,R^,Cij) + 

i=l 3 = 1 
M N 

+ Kj{{C,j + l)Ps.-i.fl,-i,c.,+i - Cy-Ps.,fl„c.,) + 

i=i j=i 

M N 

+ "Im7^J((C^J + l)Ps..i?„C.,+l - ajPs,,R^,C,,) + 
i=l i = l 

M N 

,Rj,Cij + l CijPSi,Rj,Cij) J 

^=l j=l 

C Gaussian Approximation 

We work here in some details the explicit computation for the Gaussian approximation in the specific 
case of 2 microRNAs (S'i,S'2) and 2 ceRNAs (i?i,i?2). Denoting with := a'+"LP|^,q=(i,...,i), at 
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steady state the system of equation reads: 



8 F 



d F 



d\F = 



d; . F 



F 

21,91 



F 

2:1,92 



d\F 



92 F 

Z2,9l 



F 

22,92 



F 

^91,92-' 



K-ag{dl,,F + dl^^^F) 
9s 

k..-a9{dl^,F + dl^^F) 
9s 

K-9{dl,,F + dl^^F) 
9r 

K.-9{dl,,,F + dl,q.F) 
9R 

- a9{d%,F + dl^^^^^F) 
9s 

k,,d,,F + k,,d,,F - 2ag{dl^,^^^F + d^^^^^^F) 
29s 

k„,d^,F + k^,da,F - g{d\ F + a? , „ F + ad^ ^F + adl „ „ F) 

9 + 9R+ 9S 

kq,d,,F + kAF - 9{dl.,,,F + + adj^.^^.^F + d^ j^.F) 

9 + 9R + 9S 
kAF-a9{d%,F + dl.^^^^F) 

9s 

k,,d,,F + k,,d,,F - 9{dl,,,,,,F + dl,^^^F + ad^^^^^F + adj^^^^^^F) 

9 + 9R+ 9s 

kq^d^^F + k^^dg^F - g{dl^,,,g^F + d^.^^F + ad^^^^.F + adl,,,,,,F) 



9 + 9R + 9S 



ka,d„,F - g(d'^ ^F + ^F) 

2 T? _ 21, 9^ 22, 9{ ' 



9r 

k,,dg,F + k,,d,,F - 29{dl^^^^^g^F + a3^,,„,,F) 



9 + 9R + 9S 
kq,d„^F-g(d^ ^F + d^ 2F) 

d\F = ^-^^ (16) 

9r 

Recalling that within the Gaussian approximation the partial derivatives of the third order can be 
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expressed in terms of that of lower order: 

{dl,F + d^,F)d,,F + 2d^,Fdl^^F- 
(dl.F + d.,F)d,,F + 2d.,Fdl^^^F - 



2{d,,Ffd,,F- 
2{d,,,Ffd,,F- 



^2 >92 



d-^ F 



d-^ F 

Zi,Z2,q2 

F 
F 

zi,gi,g2 



2id,,Frd,,F^ 

[dl.F + d,,F)d,,F + 2d,,Fdl^,^F ~ 2{d,,Ffd,,F 
[dl^F + d,,F)d,,F + 2d,,Fdl^^F- 2{d,,Ffd,,F^ 
{d^^F + d,2F)d,,F + 2d,,Fdl^^^F - 2id,,F)^d.,F- 

~2{d.,Ffd.,F- 
~2{d,,,Ffd,,F^ 



(d%F + d,,F)d,,F + 2d..,Fdl^^F 



F 

F 



{dt.F + d.,F)d,,F + 2d.,Fdi^^^F 
(d^fF + d,,F)d.2F + 2d,,Fdl^^^F- 2{d,,Ffd.,F 
(<9|^ + dq2F)d.2F + 2d,,Fdl^^^F ~ 2{d,,Ffd.,F 
dl ,.2 F + dl a., F + dl 9., F - 25., Fd,, F 
dl ,.2 F + af, a., F + dl 9., F - 29., F 
9.',,,, 9,2^ + dl.AF + <,,2^-2^ - 29.,F9,,F9,,F 
9f , F + dl 9,, F + dl 9., F - 29., F9,, F9,, F 



F 

Zl.92 



9^ F 



9' F 

zi,9l 



9^ F 

Zl,92 



9^ F 



9^ F 

zi,g2 



9' F 

Z2,?l 



dt „ F 



(17) 



Inserting relations (17 1 into (161 we obtain a closed system of 14 in 14 miknowns. In the general case 
of a network of N ceRNAs interacting through M miRNAs we would have a complete system of 2(A^ + 
M) + (^+*^) equations. 



D Linear noise approximation 

We use the linear noise approximation |30j in order to obtain the steady state fluctuation covariance 
matrix directly from the macroscopic system. For a general system of M miRNAs interacting with 
N mRNAs and R elementary reactions, we assign to each reaction r a propensity defined from the 
probability Qfr^ip, ^)St that a reaction r occurs in the homogeneous system volume D, in the time interval 
St. ip is the concentration vector of the M + N chemical components of the system. In the macroscopic 
limit (ri — >■ oo) the system dynamics is described by the following M + N ordinary differential equations, 

= ^iyrpfri-(pl,---,i'M+N) , (18) 

r 

where i>rp is the rp— th element of the stoichiometry matrix, i.e. it indicates the number of molecules 
by which a component p changes when an elementary reaction of type r occurs. 

F or s mall enough deviations = [Sipi, 5i}}2, SiI^m+n] from its steady state, the dynamics of Equa- 



tion ( 18 ) can be approximated by a system of linear differential equations, according to ^ = AStf), where 
A is the Jacobian matrix with elements 

= ^^rp{jr^)M>ss ■ (19) 
r=l "'^'i 

The master equation for the probability of having X = [Xi, A"2, ...,Xp, Xm+n] molecules in the 
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system at time t is then 



dP 

Itt 



{X,t) 



R N 
r— 1 p—1 



i)fr{xn-\n)p{x,t) 



(20) 



with E being a step operator with property E^'' g{..., X.^ 



, Xr, 



To obtain the hnear noise approximation 
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we expand the master equation to second order in 
after substituting each p-th component of X with Xp = U.ijjp + fi^/^Xp. Xp is the p-th component of a new 
random vector x such that the Xp is thus described as a macroscopic term il^p phis a stochastic term 
fi^/^Xp. We thus obtain a Hnear Fokkcr-Planck equation for the joint probability distribution Ii{x,t) of 



~dt 



{x,t) 



dxpli 

P,9 



IE': 



(21) 



The matrix elements apq are given by the Jacobian matrix A, while the elements bpq of the diffusion 



matrix B are defined as in 37 , 



R 



(22) 



Generally A and B may depend on time, but here we will restrict our analysis to the steady state 



case. In this way, the stationary solution of Equation (21| is the normal distribution iV(0,S). S, which 
is the covariance matrix with elements S,rpi is the solution of the matrix Lyapunov equation: 



AE + EA^ + B = 



(23) 



The covariance matrix C for the deviations in molecule number (SXi) is related to S via the relation 
C = f2S. Thus, in the linear noise approach the expected value (Xr) is approximated by Vlipr and the 
true covariance tr^p by Crp- Then, the expressions for Pearson's correlation coefficients (pXr.Xp), Fano 
factors (fx) and coefficients of variation {CVx) can be easily derived: 



PX^,Xp 
CVx^ 

fx. 



2 




(Jrr^pp 


yj Crr^pp 






{Xr) 




2 


^rr ^rr 


{Xr) 





(24) 



Therefore, thanks to Equation (23 1 the matrix C (and thus the stochastic properties of a system) can 
be directly evaluated from macroscopic parameters. 

Let's now discuss in details the specific case with two ceRNAs in interaction with one miRNAs. In 
such a system, the propensity vector / assumes the following form: 
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and the stoichiometry matrix v is given by: 



/ 1 

-1 
—a 
—a 








V 








1 

-1 

—a 
—a 











\ 



-1 





-1 




1 

-1 / 



(26) 



Thus, the Jacobian and diffusion matrices {A and B respectiveiy) foilow, 



A 



a{giiRi + 


-511^1 
— 512-R2 

/ fcsi +gsiSi +a^SiA 


a^iii^iSi 
agi2R2Si 







-9S2 -a{g2iRi +322-R2) 
-321-Ri 

— 322^^2 

agiiRiSi 
ks2 + ffS2'S'2 + o?S2B ag2iRiS2 
ag2iRiS2 ' ■ " ■ 

(722 ^^2^2 



fcj5, +RiC 




—agnSi 
—ag2iS2 
—QR — giiSi — 321S2 


aguR2Si \ 

0(722 ^?2 52 



fcjjs + R2D ) 



—agi2Si \ 
— 03225*2 


-gn — 312S1 — 5225*2 / 

(27) 



with A = giii?i +gi2i?2, B = 521^1 +522-R2, C = gn-, +giiSi+g2iS2 and D = gji^ + gi2Si+ g22S2- The 
covariance matrix elements Crp can be evaluted accordingly. In Figure |8] we plot the Pearson correlation 
coefficient of such system as a function of ceRNAl transcription rate. As it is possible to notice, Gaussian 
approximation performs better than Linear Noise approximation |311. 
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Figure 8. Comparison between Linear Noise and Gaussian approximations. (Left panel) Linear noise 
approximation, (Right panel) Gaussian approximation. Lines are analytical approximations of the 
Pearson correlation coefficient. Dots are the results of 10^ Gillespie's simulations. 
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E Figure parameters (main text) 
Figure 2 

miRNAs transcription rates: ksi = 0.05s~^ and = 0.045s~^ ; 

ceRNA2 transcription rate: ka^ = 0.155s~^ ; 

miRNA degradation rates: gs^ = gs^ = 0.0002s~^ ; 

ceRNAs degradation rates: gn^ = gn^ = 0. 0004s ~^ ; 

ceRNA-miRNA association rates: gw = gi2 = g2i = 522 = 0.0005s~^ ; 

catalyticity parameter: a = 0.1 . 

ceRNAl transcription rate is the control parameter and ranges from to 1.4s*~^ . 
Figure 3 

miRNA2 transcription rate: fcg^ = O.OSs^-'^ ; 

ceRNAs transcription rates: k^^ = 0.355s~^ and k^^ = 0.155s~^ ; 
miRNA degradation rates: gs^ = gs^ = 0.0002s~^ ; 
ceRNAs degradation rates: g^i — gn^ — 0.0004s~-'^ ; 
ceRNA-miRNA association rates: gw = g\2 = 521 = 322 = 0.0005s~^ ; 
catalyticity parameter: a = 0.1 . 

miRNAl transcription rate is the control parameter and ranges from to 0.1s*~^ . 

Figure 4 

Panel (B): 

miRNAl-10 transcription rates: kg = 0.075 + 0.01rand()s~^ ; 

ceRNA2-10 transcription rate: kji = 0.15 + 0.01rand()s~"^ ; 
miRNAl-10 degradation rates: gs = 0. 0004s ; 
ceRNAsl-10 degradation rates: gn = 0.0004s~^ ; 
miRNA-ceRNA association rates: g = 0.0006 
catalyticity parameter: a = 0.5 . 

ceRNAl transcription rate is the control parameter and every 35 hours takes the following values: 0.15s~^, 
0.5s*-\ 0.9s^-\ 0.15s^-i . 

Panel (C): 

miRNAl-9 transcription rate: ks = 0.02s~i + 0.01rand()s~i ; 
ceRNAl-10 transcription rates: kn = 0.15 + 0.01rand()s~^ ; 
miRNAl-10 degradation rates: gs = 0.0002s~^ ; 
ceRNAl-lOs degradation rates: gji = 0.0004s~-^ ; 
miRNA-ceRNA association rates: g = 0.0006s~^ ; 
catalyticity parameter: a = 0.5 . 

miRNAl transcription rate is the control parameter and every 35 hours takes the following values: 
0.02S-1, 0.1s«-i, O.Ss**-!, 0.02s«-i . 



Figure 5 

miRNAs transcription rates: ks = 0.2s~^ ; 
ceRNA2 transcription rate: ka^ = 0.155s~^ ; 
miRNA degradation rates: gsi = gs2 = 0.0003s~^ ; 
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ceRNAs degradation rates: g^^ = gji^ = 0.0004s~^ ; 
complex association rates: — k2 = 0.0005s~^ ; 
complex dissociation rates: = k2 = 0.0003s~^ ; 
complex degradation rates: 71 = 72 = O.OOOSls"^ ; 
catalyticity parameter: a = 0.1 . 

ceRNAl transcription rate is the control parameter and ranges from to 1.2s^~^ . 
Figure 6 

ceRNAl transcription rates: fc^j^ = 0.155s~^ ; 

ceRNA2o_FF-i-OA' transcription rate jumps from kj^^ = to kj^^ = 0.345 ; 

ceRNA2oN^OFF transcription rate jumps from fc^^ = 0.345 to kji^ = ; 

miRNA degradation rates: gs^ = gs-2 = 0.0002s~^ ; 

ceRNAs degradation rates: gR^ = 5_R2 — 0.0004s^^ ; 

ceRNA-miRNA association rates: gu — gi2 ~ (721 — 522 — 0.0005s~^ ; 

catalyticity parameter: a = 0.1 . 

miRNAl transcription rate is the control parameter and ranges from to O.Ss''^^ . All the other ceRNAs 
have transcription rates kji — 0.1 and all the other rates equal to ceRNAl ones. 

F Response time and experimentally testable trend 
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Figure 9. Response times and experimentally accessible parameter region. We show together the Tqn 
and Toff response times as present in Figure 6 in the main text for the case with 3 ceRNAs. The 
highlighted region corresponds to the experimentally accessible one. Increasing miRNA concentration, 
switch-off response times show a decreasing trend while switch-on a U-shaped one. The parameter 
setting is the same of Figure 6 (main text) . 
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